getwd()

brazil_vaccine <- 
  read.csv("./_3_data/brazil_vaccine.csv") %>%
  slice(-1) %>%
  left_join(., brazil_data %>% select(town_code, urb, pop) %>% mutate(town_code = as.character(town_code))) %>%
  filter(urb == 1) %>%
  select(town_code, pop, vaccine_rate, country)

indonesia_vaccine <- 
  read.csv("./_3_data/indonesia_vaccine.csv") %>%
  mutate(town_code = as.character(town_code))

vaccine_rates <-
  bind_rows(indonesia_vaccine, brazil_vaccine) %>%
  ggplot(aes(x = log10(pop), y = vaccine_rate)) +
  geom_point(color="darkgrey", shape = 21, alpha = 0.5) + 
  geom_smooth(method = "lm", color = "black") +
  stat_summary_bin(fun.y='mean', bins=20,
                   shape = 21, fill = "lightgrey",
                   color='black', size=2.5, geom='point') +
  facet_wrap(country~., scales = "free_x") +
  theme_bw() +
  theme(panel.grid.minor = element_blank(), 
        panel.grid.major.x = element_blank(),
        axis.line.y.left = element_blank(),
        legend.position = "bottom",
        strip.background = element_blank(),
        legend.title = element_blank(),
        axis.line = element_line(colour = "black"),
        panel.border = element_blank()) +
  scale_colour_grey() +
  xlab("Population, log scale") +
  ylab("Polio Vaccination Rates (%)")

ggsave("./_4_outputs/figure_oa_v.pdf", plot = vaccine_rates, width = 8, height = 5)

